#ifndef PATCH_TURBO
#define PATCH_TURBO

//FOr tesst commit

#include "grid.h"
#include <vector>
#include <string>
#include <sstream>

//Vector of variables
typedef std::vector<double> var_vector;

//Class that describes boundarys condition type
class Patch {
public:
	std::string Name;	//Name of boundary condition type
	std::vector<int> BCType; // Type of boundary condition (0 - value, 1 - flux)
	var_vector BCValue; // Value for field
	var_vector BCGrad; // Flux for field
	virtual int ParseParams (std::istringstream& ss) = 0; //Parse parameter patch string
	virtual Patch* Copy() = 0;	
	virtual var_vector Apply(Face& f, var_vector& UL) = 0;
	virtual var_vector ApplyVal(Face& f, var_vector& UL) = 0;
	//TO DO f->f.norm, delete argument Cell c
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell& c, var_vector& UL, std::vector<Vector>& grad) = 0;
};

//Sticky wall
class PatchNoSlip : public virtual Patch
{
public:
	PatchNoSlip() {
		Name = "NoSlip";
	};
	virtual int ParseParams (std::istringstream& ss) {
		return 0;
	};
	virtual Patch* Copy() {
		return new PatchNoSlip();
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[1] = -UR[1];
		UR[2] = -UR[2];
		UR[3] = -UR[3];
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[1] = 0;
		UR[2] = 0;
		UR[3] = 0;		
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell& c, var_vector& UL, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		Vector norm = f.FaceNormal/f.FaceNormal.mod();	//unity normal vector
		res[0] = grad[0] - ((norm*grad[0])*norm);	//grad_bound_cell - grad_bound_cell*n
		res[4] = grad[4] - ((norm*grad[4])*norm);
		res[1] = (norm*grad[1])*norm;				//grad_bound_cell*n for velocity
		res[2] = (norm*grad[2])*norm;
		res[3] = (norm*grad[3])*norm;
		return res;
	};
};

//Steady state
class PatchSteadyFlow : public virtual Patch
{
public:
	PatchSteadyFlow() {
		Name = "SteadyFlow";
	};
	virtual int ParseParams (std::istringstream& ss) {
		return 0;
	};
	virtual Patch* Copy() {
		return new PatchSteadyFlow();
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;		
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell& c, var_vector& UL, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		Vector norm = f.FaceNormal/f.FaceNormal.mod();	//unity normal vector
		res[0] = grad[0] - ((norm*grad[0])*norm);	//grad_bound_cell - grad_bound_cell*n
		res[1] = grad[4] - ((norm*grad[1])*norm);
		res[2] = grad[2] - ((norm*grad[2])*norm);	
		res[3] = grad[3] - ((norm*grad[3])*norm);
		res[4] = grad[4] - ((norm*grad[4])*norm);
		return res;
	};
};

//No leak slippy wall
class PatchNoLeak : public virtual Patch
{
public:
	PatchNoLeak() {
		Name = "NoLeak";
	};
	virtual int ParseParams (std::istringstream& ss) {
		return 0;
	};
	virtual Patch* Copy() {
		return new PatchNoLeak();
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		if (f.FaceSquare == 0) return UR;
		Vector VL(UL[1], UL[2], UL[3]);
		Vector VR;
		VR = VL - ( 2.0/(f.FaceNormal.mod()*f.FaceNormal.mod()) )*(VL*f.FaceNormal)*f.FaceNormal;
		UR[1] = VR.x;
		UR[2] = VR.y;
		UR[3] = VR.z;
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		Vector VL(UL[1], UL[2], UL[3]);
		Vector VR = VL - ( 1.0/(f.FaceNormal.mod()*f.FaceNormal.mod()) )*(VL*f.FaceNormal)*f.FaceNormal;
		UR[1] = VR.x;
		UR[2] = VR.y;
		UR[3] = VR.z;		
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell& c, var_vector& UL, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		Vector norm = f.FaceNormal/f.FaceNormal.mod();	//unity normal vector
		res[0] = grad[0] - ((norm*grad[0])*norm);	//grad_bound_cell - grad_bound_cell*n
		res[4] = grad[4] - ((norm*grad[4])*norm);
		Vector A = norm.x*grad[1] + norm.y*grad[2] + norm.z*grad[3]; // grad(VL*Norm)
		res[1] = grad[1] - (norm.x*A);
		res[2] = grad[2] - (norm.y*A);
		res[3] = grad[3] - (norm.z*A);
		return res;
	};
};

//Rotate velocity vector
class PatchRotate : public virtual Patch			//TO DO
{
	double alpha;
public:
	PatchRotate() {
		Name = "Rotate";
	};
	virtual int ParseParams (std::istringstream& ss) {
		ss>>alpha;
		return 0;
	};
	virtual Patch* Copy() {		
		return new PatchRotate();
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		Vector VL(UL[1], UL[2], UL[3]);
		Vector VR;
		VR.y = VL.z * sin(alpha) + VL.y * cos(alpha);
		VR.z = VL.z * cos(alpha) - VL.y * sin(alpha);
		UR[1] = VR.x;
		UR[2] = VR.y;
		UR[3] = VR.z;
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[1] = 0;
		UR[2] = 0;
		UR[3] = 0;		
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res = grad;
		res[0].y = 0;
		res[4].y = 0;
		res[1] = (-1)*U[1]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		res[2] = (-1)*U[2]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		res[3] = (-1)*U[3]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		return res;
	};
};

//Outflow P = const
class PatchOutFlowP : public Patch
{
public:
	double P;	
	double gamma;
	PatchOutFlowP(double _gamma) {
		Name = "OutFlowP";
		gamma = _gamma;
	};
	virtual int ParseParams (std::istringstream& ss) //Parse parameter patch string
	{
		ss>>P;
		return 0;
	};
	virtual Patch* Copy() {
		PatchOutFlowP* cpy = new PatchOutFlowP(gamma);
		cpy->P = P;		
		return cpy;
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		Vector v(UL[1], UL[2], UL[3]);
		v = (1.0/UL[0]) * v;			
		double cd = sqrt(gamma * (gamma - 1.0)) * sqrt(UL[4]/UL[0] - v*v/2.0);
		//double R_ = v*n - (2*cd) / (gamma - 1.0);
		if (cd < v.mod()) return UR;	//Supersonic case
		UR[4] = P/(gamma - 1.0) + (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]);
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		Vector v(UL[1], UL[2], UL[3]);
		v = (1.0/UL[0]) * v;			
		double cd = sqrt(gamma * (gamma - 1.0)) * sqrt(UL[4]/UL[0] - v*v/2.0);
		if (cd < v.mod()) return UR;
		UR[4] = 0.5*(UL[4] + P/(gamma - 1.0) + (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]));
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell& c, var_vector& UL, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		Vector norm = f.FaceNormal/f.FaceNormal.mod();	//unity normal vector
		res[0] = grad[0] - ((norm*grad[0])*norm);	//grad_bound_cell - grad_bound_cell*n   
		res[1] = grad[4] - ((norm*grad[1])*norm);
		res[2] = grad[2] - ((norm*grad[2])*norm);	//	UR=UL
		res[3] = grad[3] - ((norm*grad[3])*norm);
		Vector v(UL[1], UL[2], UL[3]);
		v = (1.0/UL[0]) * v;			
		double cd = sqrt(gamma * (gamma - 1.0)) * sqrt(UL[4]/UL[0] - v*v/2.0);
		if (cd < v.mod()){
			res[4] = grad[3] - ((norm*grad[3])*norm);		// UR=UL
		}else{

		};
		return res;
	};
};

//Outflow P = const
class PatchInFlow : public Patch			//TO DO
{
public:	
	double gamma;
	double R;	//Specific gas constant
	double P, T, vx, vy, vz; //Inflow params
	PatchInFlow(double _gamma, double _R) {
		Name = "InFlow";
		gamma = _gamma;
		R = _R;
	};
	virtual int ParseParams (std::istringstream& ss) //Parse parameter patch string
	{		
		ss>>vx>>vy>>vz>>P>>T;		
		return 0;
	};
	virtual Patch* Copy() {
		PatchInFlow* cpy = new PatchInFlow(gamma, R);	
		return cpy;
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		Vector vd(UL[1]/UL[0], UL[2]/UL[0], UL[3]/UL[0]);
		Vector n = f.FaceNormal / f.FaceNormal.mod();
		double Pd = (gamma - 1.0) * (UL[4] - (UL[1]*UL[1]+UL[2]*UL[2]+UL[3]*UL[3])/(2.0*UL[0]));
		double rod = UL[0];
		double cd2 = (gamma * Pd / rod);
		double cd = sqrt(cd2);
		double c02 = cd2 + (gamma - 1) * vd.mod() * vd.mod()/2.0;
		double R_ = vd * n - (2 * cd) / (gamma - 1);
		double cost;
		if (vd.mod() != 0) {
			cost = -vd*n/vd.mod();
		} else {
			cost = 0;
		};
		double cb = -R_ * (gamma - 1) / ((gamma - 1) * cost * cost + 2);
		//Supersonic
		if (cd2 < vd * vd) cb = sqrt(gamma * R * T);
		double tmp = ((gamma - 1) * cost * cost + 2) * c02;
		tmp /= ((gamma - 1) * R_ * R_);
		tmp -= (gamma - 1) / 2.0;
		cb *= 1 + cost * sqrt(tmp);
		double M = vd.mod() / cd;		
		double Tb = T * (cb * cb / c02);
		double Pb = P * pow((Tb/T) , gamma / (gamma - 1));
		double rob = Pb / (R*Tb);		
		double cp = R * gamma / (gamma - 1);
		double vb = sqrt(2 * cp * (T - Tb));
		//double Ptotal = Pd + (UL[1]*UL[1]+UL[2]*UL[2]+UL[3]*UL[3])/(2.0*UL[0]);		
		//double ro = Pd/(R*T);
		//double P0 = P;
		//double dP = 4.1e5;
		//double k = (vz / dP);
		Vector vr(vx, vy, vz);
		//if (M<1) {
			vr = vr * (vb / vr.mod());
		//} else {
		//	rob = P / (R * T);
		//	Pb = P;
		//};
		//vr.z = alpha * (P-Pd) * D * D;
		//vr.z = k * (P0 - Ptotal);
		//if (vr.z < 0) vr.z = 0;		
		UR[0] = rob;
		UR[1] = vr.x*rob;
		UR[2] = vr.y*rob;
		UR[3] = vr.z*rob;
		UR[4] = Pb/(gamma-1.0) + rob*(vr.x*vr.x+vr.y*vr.y+vr.z*vr.z)/2.0;
		/*Vector v(UL[1], UL[2], UL[3]);
		Vector n = (1.0/f.FaceSquare)*f.FaceNormal;
		v = (1.0/UL[0]) * v;			
		double cd = sqrt(gamma * (gamma - 1.0)) * sqrt(UL[4]/UL[0] - v*v/2.0);
		double R_ = v*n - (2*cd) / (gamma - 1.0);
		if (R_ >= 0) return UR;	//Supersonic case
		double vn = R_ + 2*sqrt(gamma*R*T) / (gamma - 1.0);*/
		//double Pd = (gamma - 1.0) * (UL[4] - (UL[1]*UL[1]+UL[2]*UL[2]+UL[3]*UL[3])/(2.0*UL[0]));
		//double Ptotal = Pd + (UL[1]*UL[1]+UL[2]*UL[2]+UL[3]*UL[3])/(2.0*UL[0]);
		/*double L = 1e-2;
		double ny = 17.2e-6;
		double D = 3.1e-3;
		double h = 40e-3;
		double alpha = 1.0 / (16*L*ny);*/
		/*
		double ro = Pd/(R*T);
		double dP = 4.1e5;
		double k = (20.0 / dP);
		Vector vr(vx,vy,vz);
		//vr.z = alpha * (P-Pd) * D * D;
		vr.z = k * (P - Ptotal);
		vr.x = 0;
		vr.y = 0;
		//vr = -vn*vr*(1.0/(vr*n));
		UR[0] = ro;
		UR[1] = vr.x*ro;
		UR[2] = vr.y*ro;
		UR[3] = vr.z*ro;
		UR[4] = Pd/(gamma-1.0) + ro*(vr.x*vr.x+vr.y*vr.y+vr.z*vr.z)/2.0;		*/
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[4] = P/(gamma - 1.0) + (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]);
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		return res;
	};
};

//Constant flow
class PatchConstant : public Patch
{
public:
	double R;
	double gamma;
	double vx,vy,vz,P,T;
	//PatchConstant(double _gamma, double _R, double MolarMass) {
	PatchConstant(double _gamma, double _R) {
		R = _R;
		gamma = _gamma;
		Name = "ConstantFlow";
	};
	virtual int ParseParams (std::istringstream& ss) //Parse parameter patch string
	{		
		ss>>vx>>vy>>vz>>P>>T;		
		return 0;
	};
	virtual Patch* Copy() {
		PatchConstant* cpy = new PatchConstant(gamma, R);	
		return cpy;
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
			double ro = P/(R*T);
			//TO DO// I think that [E] = [m^2/s^2]
			//double E = P/(gamma - 1.0) + ro * (vx*vx+vy*vy+vz*vz)/2.0;	
			//so
			double ro_E = P/(gamma - 1.0) + ro * (vx*vx+vy*vy+vz*vz)/2.0;	
			var_vector UR = UL;
			UR[0] = ro;
			UR[1] = ro*vx;
			UR[2] = ro*vy;
			UR[3] = ro*vz;
			UR[4] = ro_E;
			return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
			var_vector UR = UL;
			return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res(U.size());
		for(int i=0;i<res.size();i++) res[i] = 0.5 * grad[i];
		return res;
	};
};
//Velocity temperature constant flow
class PatchConstantP : public Patch
{
public:
	double R;
	double gamma;
	double vx,vy,vz,T;
	
	PatchConstantP(double _gamma, double _R) {
		R = _R;
		gamma = _gamma;
		Name = "ConstantFlowP";
	};
	virtual int ParseParams (std::istringstream& ss) //Parse parameter patch string
	{		
		ss>>vx>>vy>>vz>>T;		
		return 0;
	};
	virtual Patch* Copy() {
		PatchConstantP* cpy = new PatchConstantP(gamma, R);	
		return cpy;
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
			double PL = (gamma - 1.0)*(UL[4] - (UL[1]*UL[1] + UL[2]*UL[2] + UL[3]*UL[3])/(2.0*UL[0])); //Pressure inside
			double PR = PL;				
			double ro = PR/(R*T);
			double ro_E = PR/(gamma - 1) + ro*(vx*vx + vy*vy + vz*vz)/2.0;
			var_vector UR = UL;
			UR[0] = ro;
			UR[1] = ro*vx;
			UR[2] = ro*vy;
			UR[3] = ro*vz;
			UR[4] = ro_E;
			return UR;
	};//TO DO
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {			
			var_vector UR = UL;
			return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res(U.size());
		for(int i=0;i<res.size();i++) res[i] = 0.5 * grad[i];
		return res;
	};
};
/*
//Sticky wall
class PatchStickyWall : public virtual Patch
{
public:
	PatchStickyWall() {
		Name = "StickyWall";
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[1] = -UR[1];
		UR[2] = -UR[2];
		UR[3] = -UR[3];
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[1] = 0;
		UR[2] = 0;
		UR[3] = 0;		
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res = grad;
		res[0].y = 0;
		res[4].y = 0;
		res[1] = (-1)*U[1]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		res[2] = (-1)*U[2]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		res[3] = (-1)*U[3]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		return res;
	};
};

//Free border
class PatchFreeBorder : public Patch
{
public:
	PatchFreeBorder() {
		Name = "FreeBorder";
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		return res;
	};
};

//Free border P = const
class PatchFreeBorderP : public Patch
{
public:
	double P;
	double gamma;
	PatchFreeBorderP(double _P, double _gamma) {
		Name = "FreeBorderP";
		P = _P;
		gamma = _gamma;
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[4] = P/(gamma - 1.0) + (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]);
		return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
		var_vector UR = UL;
		UR[4] = P/(gamma - 1.0) + (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]);
		return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		std::vector<Vector> res = grad;
		return res;
	};
};


//Constant flow
class PatchConstant : public Patch
{
public:
	var_vector Value;
	PatchConstant(var_vector v) {
		Value = v;
		Name = "ConstantFlow";
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
			var_vector UR = Value;
			return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
			var_vector UR = Value;
			return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res(Value.size());
		for(int i=0;i<res.size();i++) res[i] = 0.5 * grad[i];
		return res;
	};
};

//Constant flow
class PatchConstantFlowP : public Patch
{
public:
	var_vector Value;
	PatchConstantFlowP(var_vector v) {
		Value = v;
		Name = "ConstantFlowP";
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
			var_vector UR = Value;
			UR[4] += (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]) - (UL[1]*UL[1] + UL[2]*UL[2] + UL[3]*UL[3]) / (2.0 * UL[0]);
			return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
			var_vector UR = Value;
			UR[4] += (UR[1]*UR[1] + UR[2]*UR[2] + UR[3]*UR[3]) / (2.0 * UR[0]) - (UL[1]*UL[1] + UL[2]*UL[2] + UL[3]*UL[3]) / (2.0 * UL[0]);
			return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res(Value.size());
		for(int i=0;i<res.size();i++) res[i] = 0.5 * grad[i];
		return res;
	};
};

//No leak and no friction
class PatchNoLeakX : public Patch
{
public:
	PatchNoLeakX() {
		Name = "NoLeakX";
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
			//TO DO
			var_vector UR = UL;
			UR[2] = -UL[2];
			UR[3] = -UL[3];
			return UR;
	};
	virtual var_vector ApplyVal(Face& f, var_vector& UL) {
			//TO DO
			var_vector UR = UL;
			UR[2] = 0;
			UR[3] = 0;
			return UR;
	};
	virtual std::vector<Vector> ApplyVisc(Face& f, Cell&c, var_vector& U, std::vector<Vector>& grad) {
		//TO DO
		std::vector<Vector> res = grad;
		res[0].y = 0;
		res[4].y = 0;
		res[1].y = 0;
		res[2] = (-1)*U[2]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		res[3] = (-1)*U[3]/((f.FaceCenter - c.CellCenter).mod()*(f.FaceNormal).mod())*f.FaceNormal;
		return res;
	};
};



class PatchFunction : public Patch
{
public:
	var_vector (*func)(Vector point);
	PatchFunction(var_vector (*f)(Vector point)) {
		func = f;
		Name = "FunctionConditions";
	};
	virtual var_vector Apply(Face& f, var_vector& UL) {
			var_vector UR = (*func)(f.FaceCenter);
			return UR;
	};
};
*/


#endif
